宗师级专家分享：地下水边界条件 源汇 等特殊问题的模拟与处理之混合抽水井的模拟、自流井模拟
首先向陈崇希老师致敬，以下来源于陈老师相关讲座，好多年之前有幸邀请陈老在年近80岁的高龄给我们分享了地下水边界条件、源汇、等特殊问题的模拟与处理讲座，现在每每想起还能感受到那位健硕而神采风扬，幽默风趣但绝对不失对学术专注、以及休息间还能跳上一段踢踏舞老者，我整理一些当时内容，分享给大家，由于内容比较多，我大概要分五篇内容，在近期逐一的分享给大家。
一
❖ 混合抽水井的模拟
❖1.概述
什么是混合井孔?依原来的定义,贯穿两个或两个以上含水层的井孔称为混合井孔(图1-1)。
这种传统上的定义实际上隐含着一个假定，即含水层内地下水是水平二维流动。
我们认为:三维流场中，常规的（非点状滤管的）抽水井或观测孔，不管其置于多层结构还是均质单层结构的含水系统，都属于混合井孔(图1-2缺) 。这是由于三维流场中,井孔滤管中不同深度的水头是不相等的，因此井管中的水要发生垂直流动，即使在不抽水的观测孔中也一样。这是混合井孔的本质所在(陈崇希等,1998b)。
由于自然界地下水普遍存在三维流(图1-3)，而井孔，特别是抽水井几乎都是非点状（滤管）井，因此混合抽水井普遍存在。
图1-1 混合井孔的传统定义
图1-3 自然界地下水普遍的三维流模式(陈崇希等,2005b）
开采地下水的井孔多采用混合井（除非水质不符合要求），以简化成井工艺，增加出水量和降低成本。因此，几乎没有无混合井的地区，特别是民间可以凿井的我国（包括台湾大量的是混合井）。既使现在开始控制打井，控制开采，但以往的混合抽水井转变为混合“观测孔”。
混合井在地下水流模型中如何刻画，国外尚未很好解决。目前主要采用由美国地调局推出，并在国际上广泛流传和应用的三维有限差分地下水流模型MODFLOW软件处理。
MODFLOW建议：“多层井的流量必须以某种方式人为地分配给每一单层，……把井流量按每层的导水系数大小分配，即 Qi /Qw=Ti /ΣT ” (1988、1996、2000)。
为了便于讨论，又不失其一般，我们以贯穿两个含水层的混合抽水井为例进行讨论。如此，上式可表述为
Q1 /Q2=T1 /T2
MODFLOW对混合井的这种处理方法没有给以理论上任何的分析、说明，缺乏理论依据，应用中与实际也不符(陈崇希等, 1992)。这是因为 (陈崇希等,2003a,2007,2008;黎明等,2008)：
（1）含水层导水系数对井孔流量的影响，不会如此简单。例如,混合井附近岩性(渗透系数)发生变化，甚至混合井打在岩性透镜体上，怎样影响井管的流量分配？含水层厚度发生变化又如何改变流量的分配？
图1-3 混合井管贯穿岩性透镜体上
（2）含水层的参数影响混合井流量的分配，导水系数只是其中一个因素，含水层的弹性给水度（储水系数）就不起作用了吗?
（3）井管的流量分配不仅与含水层的水文地质参数分布有关，还与外边界条件、井径、水泵吸水口的位置及其它抽水井的干扰等有关。
（4）按MODFLOW的“流量预分法” ，一个混合抽水井的Q1/Q2 始终是个常量(因为T1/T2为常量)。然而在模拟过程中，特别是预测时，该混合抽水井附近可加入或关闭 抽水井，在这种井群干扰下，原混合井的流量比还会保持常量吗？
（5）混合观测孔是混合抽水井的特殊情况（Qw =0）(图1-4)，对于两层混合的观测孔，其孔中水位(混合水位)必界于两含水层水位之间，即混合观测孔对于其中一含水层(例如2含水层)起抽水作用（Q2>0）, 对于另一含水层（1含水层）起注水作用（Q1<0）。如此，Q1/Q2<0。而两含水层的导水系数的比值肯定是正值，即T1/T2>0。如此，两个比值怎能相等？！
1-4 混合观测孔中的各分层的流量及漏斗曲线
陈老认为：实际上，当抽水井的流量Qw足够小时，不会改变Q2>0和Q1<0的情况,仍存在比值Q1/Q2为负值的情况，即出现负值(Q1/Q2)等于正值(T1 /T2)的不正当说法。
从理论上说，一种正确的方法，当蜕化为其特殊条件时，也必定是正确的。科学研究中，人们正是通过许多特殊来认识一般的；反过来，通过特殊来鉴别某理论的正确与否。
（6）按各层导水系数T的比例来预先人为划分各层的流量是缺乏理论依据的，因为这种做法要求各分层的有效井径相等和井壁处的水力坡度上下处处相等，这两个条件不可能人为控制，预先也不得而知。
基于上述分析，用导水系数的比值分配各层的流量是不妥的。实质上，这种方法不是模拟，而是“处理”，一种与机理不符的“处理”。
的确，Javandkl and Witherspoon(1969) 在多层井地下水流数值模拟研究中认为：当流动稳定时,在所研究的条件下,沿井壁的水力梯度是相等的。由此可得出,混合井各层的进水流量与其导水系数成正比。
但哪篇论文的井流基本模型是：和(实质上)的圆岛模型,它的结论不能随意地拓宽到一般条件。
陈老师在当时获得的几个混合井流解析模型表明(陈崇希等,2011; 陈崇希,2012)：对于傍河稳定井流；只有当各层分别为均质、等厚含水层,外边界条件相同，即两层抽水井至河流的距离相等（
另一处理混合抽水井的方法是将混合井两含水层之间的井筒给予一个很大的垂向渗透系数Kz（模拟井筒效应），通过均衡计算确定各层流量(Wen-Hsing Chiang,2001），垂向渗透系数可表示为(Neville等,2004)
式中
“防止模拟失真，提高仿真性”是水文地质模拟工作者的核心任务( 陈崇希,2003)，而地下水流系统中普遍存在的混合井，是当前国内外模拟失真的主要问题之一，应当引起我们的足够重视。
混合抽水井模拟在水质模型中的重要意义*:
思考：
上述问题促使我们去思考地下水动力学的井流分析中引入热传导源汇理论的适用性。传统地下水流模型都是纯渗流模型，因此以井孔的滤管井壁为其边界( 因为水进入滤管井壁后已不是渗流)，为此必需给出滤管井壁的边界条件——流量分布或水头分布。这是出现上述问题的根源。
❖混合抽水井模拟方法——“渗流-管流耦合模型”的提出
我们的思路是(1992)：为了避免难以给定的滤管井壁的边界条件，将井管内的水流与井周围的地下水流一起作为模拟对象，从而把模型的边界从滤管井壁移至井口(严格地讲是水泵进水口), 采用真实抽水流量或井中水位；滤管井壁的流量分布和水头分布是模拟的结果，而无需人为预先给定。所建的新模型称为“渗流-管流耦合模型” (1992)，其中渗流刻画地下水流，管流刻画井管中的水流。
对于“渗流-管流耦合模型”的求解，提出两种方法：
第一种是分别建立渗流模型和管流模型，然后采用迭代方法求解” (1992)；
第二种是把管流部分视为渗透系数很大的圆柱形透镜体，并对管流引入形式上服从达西定律的“等效渗透系数”，从而将渗流-管流模型视为含有圆柱形“透镜体”的新的“渗流”系统(1992) 。对于这个新“渗流”系统，关键问题是如何确定井管的“等效渗透系数”。
下面介绍第二种方法的具体推导.
从流体力学圆管中水流的水头损失H(伯努里)方程出发
式中：H 是水头损失(m)；f是摩擦系数；l是管长(m)；d是管内直径(m)；u是管内平均流速(m/d)；g是重力加速度(m/d2)。
当管流为层流时
式中：Re是雷诺数。而
式中：ν是流体运动粘度；μ是流体动力粘度；ρ是流体密度。
将（4-1-3）至（4-1-5）式代入（4-1-2）式，
式中：J是水力坡度；γ是流体重度。
将上式写成渗透流速的形式，对于管流，其空隙率n=1，则渗透流速V=nu=u，故有
将此式与达西定律
其中，KL定义为层流状态下管流的等效渗透系数，方程（4-1-8） 分别由陈崇希（1966）和J.Bear（1972）分别获得。
当管流呈紊流状态时，方程（4-1-2）可改写为(1992)
如果定义紊流状态下管流的等效渗透系数KN为
则
该式也具有达西定律的形式。
本来，紊流的运动规律有别于达西定律，而且不同流态区具有不同的形式。引入等效渗透系数后，将5个流态分区（1个层流区和4个紊流区）的运动规律统一为达西定律形式，且与地下水渗流定律一致。即
V=KeJ （14）
其中
为了求解水头分布，必须已知管流的水头损失ΔH，依（4-1-2）式
水头损失ΔH又依赖于f和u，f值又取决于Re(尼古拉池试验曲线 见图4-1-5 )，Re却又与u有关，即（4-1-4）式
诸要素间互相依赖,因此用迭代方法求解。为了保证必要的模拟精度，我们建议图4-1-6所示的模拟方法。
图4-1-5 Nikuradse试验曲线
图4-1-6 混合井模拟流程图
本来,涉及混合抽水井的地下水流问题属于线性-非线性流问题,其定解问题的表述是此较复杂的,引入等效渗透系数Ke 后,只要将传统的纯地下水渗流问题的数学描述,将其中的渗透系数K用等效渗透系数Ke取代,就可以了,表述简单明辽。
❖北海混合抽水试验场——混合抽水试验确定分层水文地质参数实例(1989-1992)
上述方法首次用于广西北海孔隙含水系统混合抽水试验场(陈崇希等,1992; 陈崇希等,1998) .
试验场为三个含水层夹两个弱透水层的含水系统.该试验场由间距为1.9m的双混合抽水井(Ⅰ和Ⅱ承压含水层)及8个观测孔(Ⅰ和Ⅱ承压含水层各3个,Ⅰ和Ⅱ承压含水层混合孔1个,潜水层1个)组成。观测孔至混合抽水井距离为6.94~37.14m, 混合观测孔至混合抽水井距离为10.96m。试验时间为1988年12月4日至22日:
(1) 4~10日为自然水位观测阶段,以掌握自然水位下降速率和潮汐效应变化,观测时间间隔为1h;
(2)11~16日为定流量抽水时间,两井流量之和为3394m3/d,延续时间4.5d;
(3)17~22日为停止抽水水位恢复阶段。
试验场为三个含水层夹两个弱透水层的含水系统，所建的为准三维流模型，并采用三角单元剖分的有限元数值方法。所有抽水井和观测孔都设置结点，而观测孔的实测水位经过天然动态的校正后,再用以拟合。(详见§4.13)
经过调试、优选，在不改变现场水文地质工程师提供的参数分区及参数值上下限的前提下进行拟合求参。最终的拟合动态曲线示于图4-1-7,拟合误差统计示于表4-1-1。表4-1-2 给出混合抽水井孔分层流量Q1、Q2的变化及流量以及导水系数之比的一览表,就此实例而言，MODFLOW流量分配法的误差达70%~85%。这么大的误差是不可接受的。
尽管“渗流-管流耦合模型”已用于多个实例,但由于它的重要性,我们仍对其做了室内物理模拟及其对应的数值模拟研究(陈崇希等,2004b;Chen等,2003*) ,以进一步检验该模型的仿真性。研究成果是肯定的。
图4-1-7 北海混合抽水试验 水头降深拟合曲线(据陈崇希等,1992;1998
表4-1-1 模拟水头与观测孔水位降深误差统计表 (据陈崇希 等,1992;1998)
表4-1-2 混合抽水井孔分层流量Q1、Q2的变化及流量以及导水系数比的一览表(据陈崇希等,1998)
黎明的博士论文
❖自流井的模拟
自流井是重要的人工凿井-水文地质现象。它普遍地存在于我国西北内陆盆地，冲洪积扇的前缘洼地;大型岩溶上升泉的周围以及平原区的深部承压含水层等区段。我国至今还存在大面积的地下水自流区（尽管由于大量开采地下水而缩小了其自流面积）。因此地下水流模型中如何刻画自流井是一个重要的水文地质问题。
然而，长期来未能见国外地下水流软件有符合机理的关于自流井的模拟方法。常规的纯地下水流模型已经解决不了这类问题。至今，对于数值模型中的自流井，仍是无可乃何地采用人为地给定其自流量。然而，自流井与抽水井不同，后者人们一般可以人为给定抽水流量，而井中水位是作为其响应；当井中水位降至井底时，人们就不能自由地抽取多少流量了。尽管如此，大多含有抽水井的数值模型，未能客观地将“井中水位作为抽水量的约束”来认识。然而，自流井就不能人为地给定自流量了，因为自流量是自流井对周围环境（气象、水文、地下水的开采等因素）的响应，自流井的流量是属于预测范畴，是输出信息，而不是输入信息。
能产生自流的含水层通常位于地表下深部，要通过较长的井管自流出井口。如此, 井管中水流的水头损失一般不能忽略不计，这就涉及含水层中的渗流与井中管流的耦合问题。为此我们采用“渗流-管流耦合模型” (陈崇希等,1992; 陈崇希等,1998) 可以方便地模拟出自流井的流量动态。对此，在自流井的井口处设置已知的水头边界：对于井口流速不大的自流井,采用井口水头等于井口标高即可；对于井口流速较大的自流井,则置井口水头等于井口标高加上流速水头，即
式中：为自流井井口第一类边界水头；为自流井井口标高；为自流井井口平均流速；为重力加速度。
在我国西北地下水资源-地质环境调查项目中(陈崇希等,2003d; 陈崇希等,2005b) ，此法已得到成功的应用(图4-2-1)。
这里的自流井模拟方法，在矿山中凡是自流减压排水的管井都可使用。这本来也是矿坑排水数值模拟的难点。
图4-2-1 疏勒河流域 自流井流量拟合曲线(陈崇希等,2005)（黑点为实测值；圆圈为模拟值）
❖水平井的模拟
水平井技术自30年代在前苏联和美国开始应用于油藏开采以来，已经在油藏工程中得到了广泛的应用。90年代以来，美国和其他一些发达国家陆续开始将水平井技术应用到地下水开采、地下水污染治理（地下污染物的原位抽吸系统）以及边坡排水治理工程等方面。辐射井可以看成是水平井的组合。我国在低渗透性的黄土层中开采地下水也取得很多经验。实践证明，水平井在薄含水层或低渗透性含水层抽取地下水的功能比垂直井具有明显的优势，因此引起了水文地质、石油、环境、水利工程师们的广泛关注。
有关地下水水平井井流的研究可以追朔到年Hantush 和 Papadopulos (1962)对辐射井流的研究（可视为一组水平井），他们用均匀通量分布的线汇刻画水平井管。此后，除个别外，研究者们基本上都采用线汇来刻画水平井管的作用，例如Goode等 ( 1987), Daviau等(1987), Cleveland（1994）[7], Falta（1995）[8], Zhan（1999,2001）等，对其适用性并未严格证明。
Tarshish(1992)提出考虑水平井管内水流的水头损失与平均流速的平方呈正比关系的理想模型，并限于稳定流态。然而，只有当井管中水流的流态为粗糙紊流时，井管中水流的水头损失才与平均流速呈平方关系，实际上，还存在1.0次方（层流区）、1.75次方（光滑紊流区）和两个过渡区。而且，任何一个水平井管的上游末端附近，必定存在层流区，而其他一些流态区（如层流与光滑紊流过渡区、光滑紊流区、光滑紊流与粗糙紊流过渡区以及粗糙紊流区）必然是依次出现，即只有出现前4个流态的前提下才有可能出现粗糙紊流态；另外，稳定流态不管对于油藏工程还是地下水问题都难以出现。因此, Tarshish 的方法离实际的应用还有很大的距离。均匀通量分布的线汇和等水头滤管壁面的刻画方法都未能考虑上述客观存在的复杂规律。
从流动机理来看, 水平井流就是横卧的混合井流。因此上文§4.1中提出的模拟含水层-井孔水流系统的“渗流-管流耦合模型”，同样可以用来刻画水平井流问题。陈崇希等(2002b, Chen等,2003b)研究了水平井的基本特征。
用“渗流-管流耦合模型”建立水平井的途径是: ①将管流的不同流态通过“等效渗透系数”统一地写成线性定律的形式; ②将水平井管视为渗透系数很大的园柱状透镜体，即将“含水层”视为是具有透水性很大的水平园柱状透镜体的非均质介质；③由此，可将内边界取在水平井的出水口处。如此,我们建立的地下水向有限井径水平井流动的数学模型为:
式中：
我们做了一个比较常见的河床下水平井井流的理想模型，以了解其动态特征。该含水层长度l 为116 m,宽度w为3777.4 m,厚度b为13.8 m。含水层是均质、各向同性,渗透系数K为1.0m/d;单位储水系数为0.00001 m-1。水平井位于含水平正中并与其长度l相同, 水平井的直径为0.05 m。
含水层除顶面为河流外,其它四个侧面和底面均为隔水边界;顶面河流保持10.0m的水头。初始水头为10.0m。水平井的出水口为已知水头边界,它以10.0 m/d的速度下降。要求模拟1.0d的地下水向水平井流动的动态。
我们采用多边形三维有限差分法模拟软件(陈崇希等,2005a,2007)求解上述问题。
在含水层长度l上均匀地设置了30个格点(格点距4.0 m),宽度w上不均匀地设置了25个格点(格点距0.2～1000.0 m), 厚度b上不均匀地设置了15个格点(格点距0.2～2.0 m),总共11250个格点。
模拟结果(Chen等,2003b)：水平井管单位长度的流量:模拟末时刻,由第1段的3.63 m2/d增加到最后段的9.50m2/d,后者是前者的2.6倍。这种条件下用等强度的线汇刻画水平井井管是不允许的。
平井管中水头降深s：模拟末时刻，第1结点的3.5 8 m，增加到第30结点的 10.0 m,后者水头降深是前者的2.8 倍。这种条件下用等水头的线汇刻画水平井井管也是不允许的。
等效渗透系数Ke：是随抽水延续的时间和沿水平井管的位置而变化，就抽水末时刻而言，最小值为1.46*106 m/d，最大值为6.17*107 m/d，后者是前者的41.8倍。(Chen等,2003b)这种条件下采用等阻力系数（层流的或粗糙紊流的）刻画水平井井管也是不允许的。
如果上述参数d=0.05m改为 d=0.1m，其它参数不变,则：
单位长度的流量：由第1段的22.68 m2/d逐渐增加到最后段(第29段)的92.22m2/d，后者是前者的4.07倍（ 原2.6倍）；
水头降深s：由第1结点的2.24 m逐渐增加到第30结点的 10.0 m ,后者水头降深是前者的4.47 倍。（原2.8 倍）
虽然混合抽水井流已在北海混合抽水试验场(陈崇希等,1992)得到初步验证,“渗流-管流耦合模型”提出后又应用于多个实例(陈崇希等,1995b,2003b,2003c,2003d,2004a;成建梅等，1998.),但是我们(陈崇希等,2004b)还通过砂槽物理模拟,以进一步检验包括水平井流在内的“渗流-管流耦合模型”模拟方法的可靠性。
试验的装置、步骤和结果以及数值模拟与物理试验曲线的拟合情况可见有关文献(陈崇希等,2004b) ,这里仅给出通常拟合较难的流量曲线(图4-3-1)。
图4-3-1 水平井出口流量曲线（圆点为实测值、实线为模拟值） (陈崇希等,2004b)
....................................后期持续更新
混合观测孔水位的模拟
抽水井水位的校正方法
井孔表皮效应（井周扰动效应）的处理
面井内抽水井水位降深的计算
非完整抽水井附近观测孔井水位的校正
泉的模拟
完整河流与地下水间补给、排泄的模拟
大气降雨、河渠、湖库入渗补给滞后性的刻画
潜水蒸发排泄的处理
初始水头条件的模拟
含水系统海边界问题
数值模型中流速和流量的计算问题
岩溶管道-裂隙-孔隙三重介质地下水流模拟
地下水开采-地面沉降的模拟
数据共享：
收集整理- 水文模型资料,仅供参考
数据领取：ERA5全球大气再分析数据和MSWEP多源融合降水数据
大数据！全国最新地级市和气象站空气质量和60年气候详细数据下载
最全的全国分省、市、县、乡镇行政区划矢量图（shp格式）免费下载
Ai尚研修地下水最新会议：地点西安理工大学学术交流中心
关于我们
“Ai尚研修”融合各领域的前沿技术、实践经验等，着力打造线上、线下深入高校的多元化辅助性技术学习平台，从课前到课后工作辅助“一站式”课程体系，创建学员的专属实践导师，让每一位学员达到简学践行的目的，真正的有所学、有所获、有所有。为助力科研提供重要的基础保障。
长按二维码关注Ai尚研修您随行的导师平台